Simulational Study on Dimensionality-Dependence of Heat Conduction 
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Heat conduction phenomena are studied theoretically us- 
ing computer simulation. The systems are crystal with non- 
linear interaction, and fluid of hard-core particles. Quasi-one- 
dimensional system of the size of L x x L y x L Z (L Z 3> L x ,L y ) 
is simulated. Heat baths are put in both end: one has higher 
temperature than the other. In the crystal case, the interac- 
tion potential V has fourth-order non-linear term in addition 
to the harmonic term, and Nose-Hoover method is used for 
the heat baths. In the fluid case, stochastic boundary condi- 
tion is charged, which works as the heat baths. Fourier-type 
heat conduction is reproduced both in crystal and fluid mod- 
els in three-dimensional system, but it is not observed in lower 
dimensional system. Autocorrelation function of heat flux is 
also observed and long-time tails of the form of ~ t~ d ^ 2 , where 
d denotes the dimensionality of the system, are confirmed. 

PACS number: 05.70.Ln 05.45.-a 05.20.Jj 



The normal heat conduction, which is described by the 
Fourier heat law 



J = -K- grad(T) 



(1) 



where J and T denote the heat flux and the temperature 
respectively, is one of the most fundamental nonequilib- 
rium phenomena. However the microscopic origin of it 
has not been clarified completely. Unlike the Debye's 
theory for the specific heat in equilibrium solids, this 
phenomena cannot be modeled by the harmonic chain 
jjj. In the harmonic chain, global internal temperature 
gradient is not formed, so that the thermal conductivity 
diverges. Disorder of mass of particles in the system does 
not give the saturated thermal conductivity in the ther- 
modynamic limit H . This behavior attributes to the lack 
of scattering process between modes in the system. On 
the other hand, some chaotic systems realize the Fourier 
heat law |^-|5|]. Thus the nonintegrability of the system 
which causes scattering between modes is one of crucial 
elements for the Fourier heat law. It is worth stressing 
that all these chaotic models and the Lorentz-gas model, 
which is used as the model of the heat conduction in 
metals or diluted gas are characterized by the na- 
ture that the total momentum is not conserved in the 
bulk. 

Recent studies for total momentum conserving sys- 
tem by Lepri et al. shed a light on another crucial as- 
pect, that is, dimensionality ||. They first study the 
thermodynamic behavior of energy transport in the one- 
dimensional Fermi-Pasta-Ulam (FPU) (3 lattice numer- 
ically and analytically investigated the dimensionality- 



dependence. In the numerical investigation for the one- 
dimensional case, they found numerically that the ther- 
mal conductivity diverges as N a where N is the number 
of oscillators and a is estimated roughly to be 1/2. They 
also calculated the power spectrum S(cu) of the global 
heat flux in the equilibrium and found its power-law di- 
vergence as a; -0 ' 37 in the low-frequency region. It cor- 
responds to a slow decay of autocorrelation function of 
globally averaged heat flux Cj(t) as i~ 063 . This power- 
low decay of the autocorrelation function implies the di- 
vergence of the thermal conductivity as n ~ ty q ~ ' 37 by 
applying the Green-Kubo formula p,|lO| 



1 



Cj{t)dt, 



(2) 



because the upper limit of the time integration is ex- 
pected to be proportional to the system size. Moreover 
they took an analogy with the hydrodynamical systems 
where the velocity autocorrelation function decays slowly 
as t~ d / 2 (d is the dimensionality of the system) when the 
density is high |pd|-[l5[|. In those systems, the autocor- 
relation function of heat flux also has the long-time-tails 
of the same exponent. From this analogy, they conjec- 
tured the logarithmic divergence in two-dimensional sys- 
tem and the converging thermal conductivity in three- 
dimensional system Following these studies, Hatano 
reported a diverging result in the one-dimensional di- 
atomic Toda lattice jl7| and recently Lippi and Livi con- 
firmed their prediction for the two-dimensional FPU-/3 
lattice and Lennard- Jones 6/12 lattice |18|| , 

Although the analogy with dense-fluids is plausible for 
one and two-dimensional system, there has been no di- 
rect study on cither three dimensional lattices or dense 
fluids system. Thus in the following we investigate the 
dimensionality dependencies of heat conduction on three- 
dimensional nonlinear lattice and dense fluid of hard-core 
particle. 

Our model for crystal is coupled three-dimensional 
FPU-/3 lattice with Nose-Hoover thermostats at both 
sides. Hamiltonian of this system is represented as fol- 
lows, 



N 2 



V,; 



1 = 1 
1 



i, j ;n.n 



(dqij) + -(dqij) (dqij = \q t 



(3) 



*o), (4) 



2 v V j> 4 - 

where N,l$,g denote the number of particles, natural 
length of the bond and the coefficient of the fourth or- 
der nonlinear term respectively. Mass of the particles are 



1 



uniformly 1. Here dynamical variables p i — (pfjPfiPf) 
and q { = (q^,qf,q^) are three-dimensional. The sum- 
mation over Vij takes only for nearest neighbor pairs. 
In this system, the total momentum is conserved in the 
bulk since potential is simply the function of distance 
between two particles. We will study two kinds of lat- 
tices; one is one-dimensional-like (1 x 1 x N z ) lattices, 
and the other is three-dimensional-like one (3 x 3 x N z ) 
where N z is taken to be from 8 to 256. In the follow- 
ing we will call the former as ID-lattice and the latter 
as 3D-lattice; although the 3D-lattice is actually quasi- 
one-dimensional, we will see this system maintains the 
three-dimensionality later. We take periodic boundary 
condition to the bonds vertical to ./V^-direction. Note 
that our model have no long-distance-correlations. Nose- 
Hoover thermostats jl9],^0| act on each particle in both 
ends of -/V^-direction to keep their temperature constant. 
One end is kept to be higher temperature (Tl) than the 
other (Tr). Therefore the equations of the motion of this 
system are 



Pi 



dHg 



(in the bulk) 



(5) 



-f^-CiP; (at the both ends) 



and thermostat variables move as follows; 

1). 



1 p 2 



Q 3k B T 



(6) 



Here T is the temperature of the heat bath(Tj, or Tr) and 
Q is the parameter of coupling between the thermal bath 
and the system. In this study we set Q = 1 in order that 
the response time of the thermostats (~ 1/^/Q) becomes 
the same order as the original time scale of the lattice. 
The value of other parameters are chosen as follows; 

k B = 1, T L = 152/3, T R = 8, g = 0.1, l = 100. 

Under this condition, system is fully chaotic |21]. The 
3D-lattice corresponds to the van-der-Waals solids since 
the natural length is enough long comparing with the dis- 
placement of particles which is usually smaller than 10, 
in other words, nonlinearity comes from the presence of 
the natural length is negligible, however it changes the 
strongest coupling combination of the displacements. We 
used the fourth-order predictor-corrector method to inte- 
grate the differential equations. Time step we use(10~ 3 ) 
for integration provides the total energy conservation 
with a relative accuracy 10~ 7 in the isolated system for 
the time corresponds to the longest simulation. 
We define the local temperature as 



T(i) =< p?/3 >, 



(7) 



where < • > means average over long time. The local 
heat flux jki is defined as the energy transfer per unit 
time from the fcth particle to the nearest particle j as 



Jki = P k ■ 



dV k 



hi 



(8) 



For the 3D-lattice, we will treat only fluxes flow along 
the TV^-direction bonds. Then we can represent the flux 
from the fcth particle to the next particle in ./V-direction 
simply as jk ■ The average of these heat flux jk over time 
and whole system is calculated as 



<3 > : 



1 N ' 
z< w^ jk > 

h 



(9) 



where N' is equal to (N - 1) for ID-lattice and (N - 9) 
for 3D-lattice. 

It was confirmed that the temperature profile takes a 
linear form and the effect of boundary with thermostats 
was negligible in both ends. From this fact we reasonably 
define the global thermal conductivity as 



< j > N z 
Tl-T r 



(10) 



by substituting (Tj, — Tr)/N z for the temperature gradi- 
ent in Eq.(|). 

Size-dependence of the thermal conductivity k is shown 
in FIG.|l|. In the ID-lattice, thermal conductivity di- 
verges as ~ AQ? showing a good agreement with the 
study by Lepri et al. ||. However, in the 3D-lattice, the 
result makes a sharp contrast with the lD-lattice(FIGj^). 
The thermal conductivity converges although the sys- 
tems are actually quite narrow especially for large N z 
and system become ID-lattice like in the infinite size 
limit {N z —too). 
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FIG. 1. Size dependence of the thermal conductivity in lat- 
tice systems is given in log scale. Boxes denote the one in 
ID-lattice (1 x 1 x N z ) and circles denote the one in 3D-lattice 
(3 x 3 x JV Z ). Error bars show 2a error regions. The dotted 
line is proportional to iVj" 38 which is obtained from a best fit 
by the power law for 5 points without the smallest system's 
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In order to check the analogy with hydrodynamical sys- 
tem, we also calculate temporal autocorrelation functions 
of average heat flux in the isolated systems. In this case 
the boundary conditions are periodic for all directions. 
These systems have no contact with the heat bath and 
the energy density is chosen to the value corresponding 
to (T L +T R )/2 = 88/3. In FIG.|, we can see that the au- 
tocorrelation function of ID-system decays as i -0 - 65 and 
the one of 3D-system decays i -15 asymptotically. Thus 
these results show the agreement with the Lepri et a/.'s 
conjecture. 




Time 

FIG. 2. The autocorrelation function of the isolated FPU 
lattices. Solid line is the chain system (1 x 1 x 64) and the 
dotted line is the one in the (3 x 3 x 64) lattice. Straight 
lines are obtained as the fitting in the long time region for 
power law. The gradient of the fitting line for the ID-system 
is about —0.65 and the one for the 3D-system is about —1.5. 

Next, to study the thermal conductivity of fluid, the 
system of hard-core particles is simulated. The parti- 
cles, with the uniform radius(er = 0.1) and mass(m = 1), 
collide elastically with each other. It is clear that this 
system also conserves the total momentum in the bulk. 
The dimension of the objective systems are two(Lj, xL z ), 
and three (X £ x L y x L z ) where L x and L y is fixed to 1 
and L z is taken to be from 2 to 20. The fixed width 
of the system(L x , L y ) is chosen from following condition; 
it must be large enough comparing with the radius of 
the particles, and on the other hand, it should be small 
to eliminate the hydrodynamic flow since we are inter- 
ested in the Fourier type heat conduction. The left and 
right sides of the system are the walls with temperature 
Tl = 6 and Tr = 2 (again Tl > Tr), and the boundary 
conditions for other directions are periodic. If a particle 
comes to the left or right end, where the stochastic heat 
bath is acting, velocity is randomly chosen according to 
the equilibrium distribution 

/( v ) = I(27rT)-^ i K|exp('-^V (11) 

being independent of incoming velocity. Here d denotes 



the dimension of the system and the Boltzmann constant 
is chosen to be 1. In the steady state density of the 
particles has a finite gradient; low at the hot end and 
high at the cold end. The density all over the system is 
chosen so that it will be high to model the dense fluid 
and lower than solidified density at the densest position; 
that is, 0.52 for the two-dimensional systems and 0.36 
for the three-dimensional systems relative to the close- 
packed density. 

Under these conditions, the equipartition of energy and 
the absence of the hydrodynamic flow are simulationally 
confirmed. We define the local temperature as 

n*> = s< £ *)/{ £ i), a*> 

z — e<Cr z <z+e z~€<r z <.z+e 

in the same way as the lattice case(Eq.(^)). Both in the 
two and three dimensional systems, temperature profiles 
have linear forms except for the area nearby boundary. 
For this reason, we use the temperature gradient obtained 
by the fitting in the bulk for the grad(T) in Eq.(|l|) to 
calculate the thermal conductivity of the system. The 
global heat flux is measured as the energy received by 
the cold bath per unit time in the steady state. 

Size-dependence of the thermal conductivity n is again 
investigated. In the two-dimensional system, k(L z ) di- 
verges as ~ L" where a is estimated less than 0.2. This 
result is consistent with the logarithmic divergence in 
lattice system p8| . On the other hand, as shown in 
FIG.p], k(L z ) saturates to a constant value in the three- 
dimensional system. 
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FIG. 3. Size dependence of thermal conductivity in 
three-dimensional hard core fluid system. k(L z ) saturates a 
constant value. Error bars show 2er error regions. 

The autocorrelation function of heat flux in the two 
and three-dimensional isolated systems is also observed. 
In this observation we chose the global density to the 
same value as the one of nonequilibrium steady state and 
apply the periodic boundary condition to all direction. 
The function shows power-law decay as ~ i -1 and ~ 



3 



i~ 1,5 asymptotically in the two and the three-dimensional 
system, respectively. 

In summary, using non-linear lattice and hard-core 
fluid we have confirmed normal Fourier-type heat con- 
duction in three-dimension. The ID-lattice shows the 
power-law divergence of the thermal conductivity as same 
as other one-dimensional chains. In the two-dimensional 
fluid system, thermal conductivity slowly diverges in 
thermodynamic limit. In the three-dimensional systems, 
the thermal conductivity converges in both lattice and 
fluid. We have also investigated the autocorrelation func- 
tion of the heat flux and confirmed the long time tails 
whose exponent is —d/2. These results mean that the 
divergence behavior of thermal conductivity can be un- 
derstood by the analogy with hydrodynamical system 
|^6|,Q. It is remarkable that the quasi-one-dimensional 
total momentum conserving system shows the normal 
Fourier-type heat conduction. 
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